3.8. Radial PML

by M. Wess, 2025

This Notebook is part of the dualcellspaces documentation for the addon package implementing the Dual Cell method in NGSolve, as well as part of the td_evp documentation on the implementation of time-domain methods for resonance problems.

Derivation of the PML formulation

We start from the time-domain Maxwell system in weak first-order \(E-H\) formulation

\[\begin{split} \partial_t \int_\Omega \tilde E\cdot \tilde e =\int_\Omega\mathrm{curl} \tilde H\cdot \tilde e,\\ \partial_t \int_\Omega \tilde H\cdot \tilde h = -\int_\Omega \mathrm{curl} \tilde E\cdot \tilde h \end{split}\]

for all test functions \(e,h\), where \(\Omega\) is some suitable domain.

Complex scaling

To apply a PML we change the domain of integration to some complexified mapped domain \(\gamma(\Omega)\). Additionally we introduce the covariantly mapped fields

\[\begin{split} \begin{align} E &= J^{-T}\tilde E,& e &= J^{-T}\tilde e,\\ H &= J^{-T}\tilde H,& h &= J^{-T}\tilde h, \end{align} \end{split}\]

where \(J\) is the Jacobian matrix \(J=D\gamma\). After applying the transformation rule we obtain using

\[\begin{split} \mathrm{curl}\tilde E = \mathrm{curl}\left(J^{-T} E\right) = \frac{1}{\det J}J^T\mathrm{curl} E,\\ \mathrm{curl}\tilde H = \mathrm{curl}\left(J^{-T} H\right) = \frac{1}{\det J}J^T\mathrm{curl} H, \end{split}\]

the transformed weak formulation

\[\begin{split} \partial_t \int_\Omega J^{-T} E \cdot J^{-T} e \det J=\int_\Omega\mathrm{curl} H e,\\ \partial_t \int_\Omega J^{-T} H \cdot J^{-T} h \det J = -\int_\Omega \mathrm{curl} E h. \end{split}\]

Now the above formulation with a suitable complex coordinate transformation on a suitably truncated domain is our desired yields our desired PML formulation. However, as usual the complex transformation is derived in frequency domain with a suitable dependency on the frequency. Here, we skip the detour to frequency domain and indicate this dependency nonchalantly by a dependency on time-derivation \(\partial_t\).

Namely we choose a radial scaling

\[ \gamma(x) = x + d(\|x\|)\frac{x}{\|x\|}\frac{1}{\partial t}, \]

where \(d\) vanishes in the interior domain \(\Omega_{int}=B_R(0)\) for some \(R>0\) and is a suitably chosen damping profile in \(\Omega_{ext}=\Omega\setminus\overline{\Omega_{int}}=\overline{B_R(0)}^c\). A simple case is e.g., given by

\[\begin{split} d(r) = \begin{cases} 0,&r<R\\ \alpha (r-R),&r\geq R \end{cases} \end{split}\]

for a given damping constant \(\alpha>0\).

Thus, we have for the scaling Jacobian \(J\)

\[ \begin{align} J(x) &= I + \frac{d(\|x\|)}{\|x\|}I\frac{1}{\partial_t}+\left(\frac{d'(\|x\|)}{\|x\|^2}-\frac{d(\|x\|)}{\|x\|^3}\right)xx^T\frac{1}{\partial_t} \end{align} \]

We may rewrite the above using the (orthogonal) projections onto \(x\) and \(x^\perp\) using

\[ \begin{align} \Pi_{x} &= \frac{ xx^T}{\|x\|^2},&\Pi_x^\perp = I-\frac{xx^T}{\|x\|^2} \end{align} \]

to obtain

\[ J(x) = \left(1+d'(\|x\|)\frac{1}{\partial_t}\right)\Pi_x+\left(1+\frac{d(\|x\|)}{\|x\|}\frac{1} {\partial_t}\right)\Pi_x^\perp \]

and thus

\[ J^{-T}=\frac{\partial_t}{\partial_t+d'}\Pi_x+\frac{\partial_t \|x\|}{\|x\|\partial_t+d}\Pi_x^\perp \]

\[\begin{split} \begin{align} \partial_t J^{-1}J^{-T}\det J &= \frac{1}{\partial_t+d'}\frac{(\|x\|\partial_t+d)^2}{\|x\|^2}\Pi_x+(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t+\frac{2d\partial_t}{\|x\|}\frac{1}{\partial_t+d'}+\frac{d^2}{\|x\|^2}\frac{1}{\partial_t+d'}-\frac{d'\partial_t}{\partial_t+d'}\right)\Pi_x+(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t +\frac{2d}{\|x\|}-\frac{2dd'}{\|x\|}\frac{1}{\partial_t+d'} +\frac{d^2}{\|x\|^2}\frac{1}{\partial_t+d'} -d'+\frac{d'd'}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t-d'+\frac{2d}{\|x\|} +\left(-\frac{2dd'}{\|x\|} +\frac{d^2}{\|x\|^2} +d'd' \right)\frac{1}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp\\ &=\left(\partial_t-d'+\frac{2d}{\|x\|} +\frac{\left(d'-\frac{d}{\|x\|} \right)^2}{\partial_t+d'} \right)\Pi_x +(\partial_t+d')\Pi_x^\perp \end{align} \end{split}\]

The full system

introducing the auxiliary unknowns

\[E_1:= \frac{d'-\frac{d}{\|x\|}}{\partial_t+d'} E\]

and similarly for \(H_1\) leads to

\[\begin{split} \partial_t E_1 = (d'-\frac{d}{\|x\|})E-d'E_1\\ \partial_t H_1 = (d'-\frac{d}{\|x\|})H-d'H_1. \end{split}\]

Thus, we obtain the full PML system

\[\begin{split} \begin{align} \partial_t \int_\Omega E\cdot e &=\int_\Omega(d'-\frac{2d}{\|x\|})\Pi_x E\cdot e+\int_\Omega(-d'+\frac{d}{\|x\|})\Pi_x E_1\cdot e-\int_\Omega d' \Pi_x^\perp E\cdot e+\int_\Omega\mathrm{curl} H e,\\ \partial_t\int_\Omega E_1 \cdot e_1 &= \int (d'-\frac{d}{\|x\|})E\cdot e_1-\int d' E_1\cdot e_1\\ \partial_t \int_\Omega H\cdot h &= \int_\Omega(d'-\frac{2d}{\|x\|})\Pi_x H\cdot h+\int_\Omega(-d'+\frac{d}{\|x\|})\Pi_x H_1\cdot h-\int_\Omega d' \Pi_x^\perp H\cdot h-\int_\Omega \mathrm{curl} E h\\ \partial_t\int_\Omega H_1\cdot h_1 &= \int (d'-\frac{d}{\|x\|})H\cdot h_1-\int d' H_1\cdot h_1 \end{align} \end{split}\]

Implementation

We start by doing the necessary imports

from ngsolve import *
import dualcellspaces as dcs
from ngsolve.webgui import Draw
from netgen.occ import *
import numpy as np
import matplotlib.pyplot as pl

R = 0.8     # inner sphere radius
R_pml = 1.1     # outer sphere radius

maxh = 0.06

order = 0

alpha = 10


# Create outer and inner spheres (both centered at origin)
outer_sphere = Sphere(Pnt(0, 0, 0), R_pml)
inner = Sphere(Pnt(0, 0, 0), R)

pml = outer_sphere-inner
pml.name = "pml"
inner.name = "inner"

geo = OCCGeometry(Glue([pml,inner]))
mesh = Mesh(geo.GenerateMesh(maxh = maxh))

#mesh.Curve(order)
Draw(mesh)
print(mesh.GetMaterials())

('pml', 'inner')

fes_E = dcs.HCurlDualCells3D(mesh, order=order)
fes_H = dcs.HCurlPrimalCells3D(mesh,order=order)

print("total DoFs:",fes_E.ndof+fes_H.ndof)
gfE, gfH = GridFunction(fes_E), GridFunction(fes_H)
gfE1, gfH1 = GridFunction(fes_E), GridFunction(fes_H)


#integral symbols with special integration rules
dxE = dx(intrules=fes_E.GetIntegrationRules())
dxH = dx(intrules=fes_H.GetIntegrationRules())

dxw = dx(intrules=dcs.GetIntegrationRules(2*order+4))
dSw = dx(element_boundary=True,intrules=dcs.GetIntegrationRules(2*order+4))


#mixed bilinear form
E,dE = fes_E.TnT()
H,dH = fes_H.TnT()

normal = specialcf.normal(3)
bf_mixed = BilinearForm(E*curl(dH)*dxw+E*Cross(dH,normal)*dSw, geom_free=True).Assemble().mat

total DoFs: 489292

pmldofs_E = fes_E.GetDofs(mesh.Materials('pml'))
pmldofs_H = fes_H.GetDofs(mesh.Materials('pml'))

# coefficients for mass matrices
vx = CF((x,y,z))
r = Norm(vx)

Pi = OuterProduct(vx/r,vx/r)
Pi_perp = Id(3)-Pi


d = mesh.MaterialCF({"inner" : 0, "pml" : R+(r-R)*alpha})

d_ = d.Diff(r)


#Draw(d,mesh,draw_surf=False, clipping={"y":0, "z":-1},
#      settings = {"Objects" : {"Clipping Plane" : True}}, euler_angles=[0,40,0])

#Draw(d_,mesh,draw_surf=False, clipping={"y":0, "z":-1},
#      settings = {"Objects" : {"Clipping Plane" : True}}, euler_angles=[0,40,0])

#prepare mass operators
with TaskManager():
    # (inverse) mass
    massH_inv = fes_H.Mass(1).Inverse()
    massE_inv = fes_E.Mass(1).Inverse()
    
    # inverse mass in pml only   
    massH_inv_pml = fes_H.Mass(1).Inverse(freedofs = pmldofs_H)
    massE_inv_pml = fes_E.Mass(1).Inverse(freedofs = pmldofs_E)    
    
    # mass matrices for pml damping
    dampEe = BilinearForm((d_-2*d/r)*(Pi*E)*dE*dxE-d_*Pi_perp*E*dE*dxE).Assemble().mat
    dampE1e = BilinearForm((-d_+d/r)*(Pi*E)*dE*dxE).Assemble().mat
    dampEe1 = BilinearForm((d_-d/r)*E*dE*dxE).Assemble().mat
    dampE1e1 = BilinearForm(-d_*E*dE*dxE).Assemble().mat

    dampHh = BilinearForm((d_-2*d/r)*(Pi*H)*dH*dxH-d_*Pi_perp*H*dH*dxH).Assemble().mat
    dampH1h = BilinearForm((-d_+d/r)*(Pi*H)*dH*dxH).Assemble().mat
    dampHh1 = BilinearForm((d_-d/r)*H*dH*dxH).Assemble().mat
    dampH1h1 = BilinearForm(-d_*H*dH*dxH).Assemble().mat
    

def estimate_tau(mat, maxsteps = 2000, tol = 1e-9):   
    vec = mat.CreateColVector()
    vec.SetRandom()
    vec*=1/vec.Norm()
    tmp = vec.CreateVector()
    lam = 0
    
    for i in range(maxsteps):
        #print(i,end='\r')
        tmp.data = mat * vec
        
        lamnew = InnerProduct(tmp,vec)
        tau = 2/sqrt(lamnew)
        n = tmp.Norm()
        res = (tmp-lam*vec).Norm()/lamnew/n
        #print(res)
        tmp *= 1/n
        if res<tol: return tau
        vec.data = tmp
        lam = lamnew
    print("did not converge, last res = ",res)
    return tau

tau = estimate_tau(massE_inv@bf_mixed.T@massH_inv@bf_mixed)
tau *= 0.9
print("estimated tau = ", tau)

estimated tau =  0.010683411757906781

#set initial data

a = exp(-r**2*20)


E0 = CF((a.Diff(y),-a.Diff(x),0))  #E0 = curl(a*(0,0,1))
rhs = LinearForm(E0*dE*dxE).Assemble().vec
gfE.vec.data = massE_inv*rhs

#Draw solution
#scenee = Draw (gfE.Operator("altshape"), mesh, "E", order=2, draw_surf=False, clipping={"y":0, "z":-1},
#               min=-0.5,max=0.5,autoscale=False,points=dcs.GetWebGuiPoints(2), vectors=True,
#               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])

# or save solution in multivector for animation
#gfE_anim = GridFunction(fes_E,multidim=0)

from time import time

drawevery = 70
tend = 1.1

t = 0.
i = 0

gfH.vec[:]=0
#gfE.vec[:]=0

gfH1.vec[:]=0
gfE1.vec[:]=0

tmpE = gfE.vec.CreateVector()
tmpH = gfH.vec.CreateVector()

tmpE_ = gfE.vec.CreateVector()
tmpH_ = gfH.vec.CreateVector()



now = time()
timepassed = 0


damp = True

with TaskManager():
    gfH.vec.data += -tau/2*massH_inv@bf_mixed*gfE.vec     #initial data in pml is zero, so first step does not need damping terms
    while t<tend:
        if i%drawevery == 0:
            timepassed += time()-now
            #scenee.Redraw() # for immediate drawing
            #gfE_anim.AddMultiDimComponent(gfE.vec) #for saving to multidim function
            Draw (gfE.Operator("altshape"), mesh, "E0", order=1, draw_surf=False, clipping={"y":0, "z":-1},
               min=-1,max=1,autoscale=False,points=dcs.GetWebGuiPoints(1), 
               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])            
            print(" time = {}, norm = {}, step = {}, {:e} dofs/s".format(t,gfE.vec.Norm(),i,i*(fes_E.ndof+fes_H.ndof)/timepassed),end="")

            now = time()
        t+=tau       
        i+=1

        if damp:
            tmpE.data = dampE1e1*gfE1.vec+dampEe1*gfE.vec
            tmpE_.data = bf_mixed.T*gfH.vec+dampEe*gfE.vec+dampE1e*gfE1.vec
            gfE.vec.data += tau*massE_inv*tmpE_
            gfE1.vec.data += tau*massE_inv_pml*tmpE


            tmpH.data = dampH1h1*gfH1.vec+dampHh1*gfH.vec
            tmpH_.data = -bf_mixed*gfE.vec+dampHh*gfH.vec+dampH1h*gfH1.vec
            gfH.vec.data += tau*massH_inv*tmpH_
            gfH1.vec.data += tau*massH_inv_pml*tmpH
        
        else:
            gfE.vec.data += tau*massE_inv@bf_mixed.T*gfH.vec
            gfH.vec.data += -tau*massH_inv@bf_mixed*gfE.vec


        

 time = 0.0, norm = 4.573211163926184, step = 0, 0.000000e+00 dofs/s

 time = 0.747838823053474, norm = 2.750803648279417, step = 70, 1.382374e+07 dofs/s

#print("finshed with an average of {:e} dofs/s".format(i*(fes_E.ndof+fes_H.ndof)/timepassed),end="")

#Draw (gfE_anim, mesh, "E0", order=2, draw_surf=False, clipping={"y":0, "z":-1},
#               min=-1,max=1,autoscale=False,points=dcs.GetWebGuiPoints(1), 
#               animate = True,
#               settings = {"eval": 0, "Objects" : {"Clipping Plane" : True}}, euler_angles=[0,30,0])